package ffte.algorithm

/**
 * FFT Generator - Core FFT algorithm generation engine
 *
 * This object provides the main API for generating FFT components with different
 * decomposition strategies including prime factorization, mixed-radix, and Winograd methods.
 *
 * Supported Algorithms:
 * - PSB (Prime Space Base): Uses prime factor decomposition for coprime factors
 * - MSB (Mixed Space Base): Uses mixed-radix decomposition for general cases
 * - WPB (Winograd Prime Base): Uses Winograd algorithms for prime factors
 *
 * Generation Strategy:
 * The generator automatically selects the optimal decomposition method based on:
 * - Input size N and its prime factorization
 * - Whether Winograd or standard DFT methods are preferred
 * - Memory and computational constraints
 *
 * Algorithm Selection Methods:
 * - Winograd: Use Winograd algorithms for prime factors when possible
 * - FlattenDFT: Use standard DFT decomposition methods
 *
 * @example
 * // Generate 64-point FFT using default method
 * val fft64 = FFTGen.gen[Flat](64)
 *
 * // Generate 30-point FFT with Winograd preference
 * FFTGen.Winograd
 * val fft30 = FFTGen.gen[Mixed](30)
 *
 * // Generate reverse FFT structure
 * val revFFT = FFTGen.rgen[Flat](128)
 *
 * @see gen for forward FFT generation
 * @see rgen for reverse FFT generation
 * @see choose_factor for factor selection algorithm
 */
object FFTGen {
    /**
     * Generic Fourier Transform trait
     *
     * Defines the basic interface for all Fourier transform implementations.
     * Provides the size parameter N that represents the transform length.
     *
     * This trait serves as the foundation for all FFT decomposition strategies
     * including DFT, PSB, MSB, and WPB implementations.
     *
     * @tparam N Transform size (number of points)
     */
    trait genericFT {
        /** Transform size - number of points in the FFT */
        def N: Int
    }

  /**
     * Discrete Fourier Transform (DFT) case class
     *
     * Represents a basic DFT implementation using the direct formula:
     * X[k] = Σ[n=0 to N-1] x[n] * exp(-2πi*k*n/N)
     *
     * This is the fundamental building block for all FFT decompositions.
     * Larger FFTs are decomposed into smaller DFT operations.
     *
     * Mathematical Foundation:
     * The DFT transforms a sequence of N complex numbers from the time domain
     * to the frequency domain. Each output X[k] represents the amplitude and
     * phase of the frequency component k.
     *
     * @param N Size of the DFT (number of points)
     *
     * @example
     * val dft8 = DFT(8)  // Create 8-point DFT
     * val dft16 = DFT(16) // Create 16-point DFT
     */
    case class DFT(val N: Int) extends genericFT {
        /**
         * String representation for debugging and logging
         * @return String in format "dft_N" where N is the transform size
         */
        override def toString = s"dft_$N"
    }

    /**
     * Combined Fourier Transform trait
     *
     * Represents FFT decompositions that combine two smaller FFT operations.
     * This is the base trait for multi-stage FFT algorithms like PSB, MSB, and WPB.
     *
     * Decomposition Strategy:
     * Large FFTs are decomposed into two smaller FFTs using the Cooley-Tukey algorithm:
     * X[k] = Σ[m=0 to M-1] Y[m] * W_N^(k*m) where Y[m] is the m-th sub-transform
     *
     * Architecture Pattern:
     * - First stage (ib): Inner transform operating on groups of data
     * - Second stage (ob): Outer transform combining results from first stage
     * - Twiddle factors: Applied between stages for proper frequency mapping
     *
     * @tparam T Implementation method type (Flat, Mixed, Streamed)
     */
    trait combineFT[T <: Method] extends genericFT {
        /** FFT decomposition structure containing stage information */
        val p: ffter[T]

        /** Inner (first) stage FFT - operates on data groups */
        val ib = p.first.fft

        /** Outer (second) stage FFT - combines first stage results */
        val ob = p.second.fft
    }

    /**
     * Prime Space Base FFT case class
     *
     * Implements FFT decomposition using the Prime Space Base (PSB) algorithm.
     * Used when the two factors are coprime (gcd = 1), enabling efficient
     * memory access patterns and reduced twiddle factor complexity.
     *
     * Algorithm Principle:
     * When N = N1 * N2 where gcd(N1, N2) = 1, the DFT can be computed using:
     * 1. N1-point DFTs on N2 groups of data
     * 2. N2-point DFTs on N1 groups of results
     * 3. Simplified twiddle factors due to coprimality
     *
     * Advantages:
     * - Optimal memory access patterns
     * - Reduced twiddle factor complexity
     * - Efficient for coprime factorizations
     *
     * @param p FFT decomposition structure with coprime factors
     * @tparam T Implementation method type
     */
    case class psbFFT[T <: Method](p: ffter[T]) extends combineFT[T] {
        /** Total transform size from decomposition structure */
        val N = p.rN

        /**
         * String representation for debugging
         * @return String showing PSB decomposition with inner and outer transforms
         */
        override def toString = s"psb f_$ib,s_$ob"
    }

    /**
     * Mixed Space Base FFT case class
     *
     * Implements FFT decomposition using the Mixed Space Base (MSB) algorithm.
     * Used for general factor combinations where factors may share common divisors.
     *
     * Algorithm Principle:
     * When N = N1 * N2 (general case), the DFT is computed using:
     * 1. N1-point DFTs on N2 groups of data
     * 2. N2-point DFTs on N1 groups of results
     * 3. Full twiddle factor multiplication
     *
     * Advantages:
     * - Works for any factorization
     * - General-purpose decomposition
     * - Handles non-coprime factors
     *
     * Tradeoffs:
     * - More complex twiddle factors than PSB
     * - Less optimal memory access patterns
     *
     * @param p FFT decomposition structure with general factors
     * @tparam T Implementation method type
     */
    case class msbFFT[T <: Method](p: ffter[T]) extends combineFT[T] {
        /** Total transform size from decomposition structure */
        val N = p.rN

        /**
         * String representation for debugging
         * @return String showing MSB decomposition with inner and outer transforms
         */
        override def toString = s"msb f_$ib,s_$ob"
    }

    /**
     * Winograd Prime Base FFT case class
     *
     * Implements FFT decomposition using Winograd's prime factor algorithms.
     * Specifically designed for prime-sized FFTs using mathematical optimizations
     * that reduce the number of multiplications.
     *
     * Algorithm Principle:
     * Winograd algorithms reduce multiplication count through:
     * 1. Pre-computation of convolution kernels
     * 2. Use of number-theoretic transforms
     * 3. Optimized mapping for prime sizes
     *
     * Mathematical Foundation:
     * Based on Winograd's work on minimal multiplication algorithms
     * for computing short convolutions and DFTs of prime length.
     *
     * Advantages:
     * - Minimal multiplication count for prime sizes
     * - Optimized for hardware implementation
     * - Reduced computational complexity
     *
     * @param N Transform size (must be prime)
     * @param p FFT decomposition structure for Winograd implementation
     * @tparam T Implementation method type
     */
    case class wpbFFT[T <: Method](val N: Int, p: ffter[T]) extends combineFT[T] {
        /**
         * String representation for debugging
         * @return String showing Winograd decomposition with inner and outer transforms
         */
        override def toString = s"wino f_$ib,s_$ob"
    }

    /**
     * FFT implementation method trait
     *
     * Defines the interface for different FFT implementation strategies.
     * Each method represents a different hardware implementation approach
     * with different resource usage, timing, and performance characteristics.
     *
     * Implementation Types:
     * - Flat: Combinatorial implementation with single-cycle computation
     * - Mixed: Mixed-radix implementation with internal pipelining
     * - Streamed: Fully pipelined streaming implementation
     *
     * Method Selection Criteria:
     * - Resource constraints (DSP usage, memory requirements)
     * - Timing requirements (latency, throughput)
     * - Power consumption considerations
     * - Clock frequency targets
     *
     * @tparam T Specific implementation method type
     */
    trait Method {
        /**
         * Apply the implementation method to generate hardware
         *
         * This is the core method that translates the abstract FFT specification
         * into concrete hardware implementation based on the selected method.
         *
         * @param entity The FFT specification to implement
         * @param prefix String prefix for generated signal names
         * @return Hardware implementation component
         */
        def apply(entity: genericFT, prefix: String): Any
    }
    
    /**
     * FFT generator trait - Core decomposition engine
     *
     * This trait represents the fundamental FFT decomposition structure.
     * It defines how a large FFT is broken down into smaller FFT operations
     * and manages the data flow between stages.
     *
     * Architecture Pattern:
     * Each ffter represents a two-stage decomposition:
     * - First stage: Processes groups of input data
     * - Second stage: Combines results from first stage
     * - Load/Store tables: Manage data reordering between stages
     *
     * Mathematical Foundation:
     * Based on the Cooley-Tukey FFT decomposition formula:
     * X[k] = Σ[n=0 to N-1] x[n] * W_N^(k*n)
     * where W_N = exp(-2πi/N) is the twiddle factor
     *
     * Data Flow Management:
     * - Load table: Defines how to read input data for first stage
     * - Store table: Defines how to write output data from second stage
     * - Index mapping: Handles the permutation of data between stages
     *
     * @tparam T Implementation method type
     */
    trait ffter[T <: Method] {
        /** First (inner) stage FFT processor */
        def first: ffter[Method]

        /** Second (outer) stage FFT processor */
        def second: ffter[Method]

        /**
         * Total transform size
         *
         * Computed as the product of both stage sizes.
         * This represents the overall FFT point count.
         */
        def rN: Int = first.rN * second.rN

        /** FFT specification for this decomposition level */
        def fft: genericFT

        /**
         * Load table - input data permutation indices
         *
         * Defines how input data should be reordered before processing
         * by the first stage. Each index specifies which input sample
         * should be fed to which input of the first stage.
         *
         * Size: rN elements
         * Values: 0 to rN-1 (input sample indices)
         */
        val load_tab: Array[Int]

        /**
         * Store table - output data permutation indices
         *
         * Defines how output data from the second stage should be reordered
         * to produce the final FFT output. Each index specifies where
         * each output should be placed in the final result.
         *
         * Size: rN elements
         * Values: 0 to rN-1 (output sample indices)
         */
        val store_tab: Array[Int]

        /**
         * Calculate load index for two-dimensional data access
         *
         * Converts 2D coordinates (f,s) to 1D load index.
         * Used for accessing input data in the correct order.
         *
         * @param f First dimension index (0 to first.rN-1)
         * @param s Second dimension index (0 to second.rN-1)
         * @return Load table index for input data
         */
        def load(f: Int, s: Int): Int = f * second.rN + s

        /**
         * Calculate store index for two-dimensional data access
         *
         * Converts 2D coordinates (f,s) to 1D store index.
         * Used for storing output data in the correct order.
         *
         * @param f First dimension index (0 to first.rN-1)
         * @param s Second dimension index (0 to second.rN-1)
         * @return Store table index for output data
         */
        def store(f: Int, s: Int): Int = f + s * second.rN

        /**
         * Generate hardware implementation
         *
         * Creates the actual hardware component for this FFT decomposition
         * using the specified implementation method.
         *
         * @param m Implementation method to use
         * @return Hardware implementation component
         */
        def gen(m: T) = m(fft, "top")

        /**
         * Reverse generation flag
         *
         * Indicates whether this structure is configured for forward
         * or reverse FFT generation. Affects twiddle factor signs
         * and data ordering.
         */
        def rgen: Boolean

        /**
         * Debug dump of permutation tables
         *
         * Creates a string representation of both load and store tables
         * for debugging and analysis purposes.
         *
         * @return String in format "load1,load2,...|store1,store2,..."
         */
        def tab_dump = load_tab.map(_.toString).reduce((a,b) => a + "," + b) + "|" + store_tab.map(_.toString).reduce((a,b) => a + "," + b)
    }

    /**
     * Base case FFT - 1-point DFT
     *
     * This object serves as the termination case for FFT decomposition.
     * A 1-point FFT is essentially a no-op operation since the DFT
     * of a single sample is the sample itself.
     *
     * Mathematical Foundation:
     * For N=1, the DFT formula reduces to: X[0] = x[0] * exp(0) = x[0]
     * This makes the 1-point FFT the simplest possible case.
     *
     * Role in Decomposition:
     * - Serves as the base case for recursive FFT decomposition
     * - Terminates the factorization tree
     * - Provides a well-defined endpoint for all algorithms
     *
     * Implementation Notes:
     * - First and second stages reference itself (recursive definition)
     * - Load and store tables are trivial (single element)
     * - No actual computation required
     */
    case object fft1 extends ffter[Method] {
        /** Recursive reference to self for first stage */
        val first = fft1

        /** Recursive reference to self for second stage */
        val second = fft1

        /**
         * Transform size for base case
         * @return 1 (single point transform)
         */
        override def rN = 1

        /**
         * FFT specification for base case
         * @return 1-point DFT specification
         */
        def fft: genericFT = DFT(1)

        /**
         * Load table for base case
         * Single element maps to index 0
         */
        val load_tab: Array[Int] = Array(0)

        /**
         * Store table for base case
         * Single element maps to index 0
         */
        val store_tab: Array[Int] = Array(0)

        /**
         * Reverse generation flag
         * Base case doesn't distinguish between forward/reverse
         */
        val rgen = false
    }

    /**
     * Prime Base FFT class
     *
     * Implements FFT for prime-sized transforms using direct DFT computation.
     * Used when the transform size N is prime and cannot be further factored.
     *
     * Algorithm Strategy:
     * For prime N, no further factorization is possible, so we use
     * the direct DFT formula: X[k] = Σ[n=0 to N-1] x[n] * exp(-2πi*k*n/N)
     *
     * Mathematical Properties:
     * - Prime numbers have no non-trivial factors
     * - Cannot apply Cooley-Tukey decomposition
     * - Direct DFT is the only approach (except Winograd)
     *
     * Implementation Characteristics:
     * - Load table: Identity mapping (0,1,2,...,N-1)
     * - Store table: Identity mapping (0,1,2,...,N-1)
     * - No data reordering required
     * - Single-stage computation
     *
     * Performance Considerations:
     * - O(N²) computational complexity
     * - Suitable for small prime sizes (typically N < 32)
     * - For larger primes, consider Winograd algorithms
     *
     * @param N Transform size (must be prime)
     * @param rgen Whether to use reverse generation
     * @tparam T Implementation method type
     *
     * @throws AssertionError when N is not prime
     *
     * @example
     * val prime5 = pb[Flat](5, false)  // 5-point prime FFT
     * val prime7 = pb[Mixed](7, true)  // 7-point prime FFT, reverse
     */
    case class pb[T <: Method](N: Int, rgen: Boolean) extends ffter[T] {
        // Validate that N is prime
        assert(Prime.isPrime(N), s"pb need $N is prime")

        /** First stage references itself (single-stage computation) */
        val first = this.asInstanceOf[ffter[Method]]

        /** Second stage is the base case (no actual second stage) */
        val second = fft1

        /**
         * Transform size
         * @return N (prime transform size)
         */
        override def rN = N

        /** FFT specification - direct DFT of size N */
        def fft = DFT(N)

        /**
         * Load index calculation
         *
         * For prime base FFT, data access is direct with no reordering.
         * The load operation simply returns the first index.
         *
         * @param p First dimension index (ignored for prime case)
         * @param s Second dimension index (ignored for prime case)
         * @return Direct index p
         */
        override def load(p: Int, s: Int): Int = p

        /**
         * Store index calculation
         *
         * For prime base FFT, data storage is direct with no reordering.
         * The store operation simply returns the first index.
         *
         * @param p First dimension index (ignored for prime case)
         * @param s Second dimension index (ignored for prime case)
         * @return Direct index p
         */
        override def store(p: Int, s: Int): Int = p

        /**
         * String representation for debugging
         * @return String showing prime FFT size
         */
        override def toString = s"prime fft $N"

        /**
         * Load table - identity permutation
         *
         * For prime FFT, input data is accessed in natural order
         * with no reordering required.
         */
        val load_tab: Array[Int] = (for (i <- 0 until N) yield i).toArray

        /**
         * Store table - identity permutation
         *
         * For prime FFT, output data is stored in natural order
         * with no reordering required.
         */
        val store_tab: Array[Int] = (for (i <- 0 until N) yield i).toArray
    }

    /**
     * Prime Space Base FFT decomposition class
     *
     * Implements PSB (Prime Space Base) FFT decomposition for coprime factors.
     * Used when N = N1 * N2 where gcd(N1, N2) = 1, enabling efficient
     * computation using number theory properties.
     *
     * Mathematical Foundation:
     * Based on the Prime Factor FFT algorithm which exploits the Chinese
     * Remainder Theorem when factors are coprime. The key insight is that
     * for coprime N1 and N2, there exists a one-to-one mapping between
     * the single index k (0 to N-1) and the pair of indices (k1, k2).
     *
     * Chinese Remainder Theorem Application:
     * For coprime N1, N2, each k ∈ [0, N-1] can be uniquely represented as:
     * k ≡ k1 (mod N1) and k ≡ k2 (mod N2)
     * This enables efficient index mapping and twiddle factor reduction.
     *
     * Data Reordering Strategy:
     * - Load table: Maps input indices to two-dimensional grid
     * - Store table: Maps two-dimensional results back to output indices
     * - Uses Chinese Remainder Theorem for index conversion
     *
     * @param first First stage FFT processor (size N1)
     * @param second Second stage FFT processor (size N2)
     * @param rgen Whether to use reverse generation mode
     * @tparam T Implementation method type
     *
     * @throws AssertionError when first.rN and second.rN are not coprime
     *
     * @example
     * // Create PSB for 15-point FFT (3 * 5, coprime)
     * val stage3 = pb[Flat](3, false)
     * val stage5 = pb[Flat](5, false)
     * val psb15 = psb[Flat](stage3, stage5, false)
     */
    case class psb[T <: Method](first: ffter[Method], second: ffter[Method], rgen: Boolean) extends ffter[T] {
        /** Total transform size N = N1 * N2 */
        val N = rN

        // Validate that factors are coprime
        assert(Prime.gcd(first.rN, second.rN) == 1, "psb factors must be coprime")

        /** FFT specification - use PSB algorithm */
        def fft = psbFFT(this)

        /** Chinese Remainder Theorem base for index mapping */
        val base = Prime.Base(first.rN, second.rN)

        /**
         * Load table generation for PSB decomposition
         *
         * Maps input indices to two-dimensional grid access pattern.
         * Two different strategies based on generation mode:
         *
         * Forward generation (rgen=false): Use Chinese Remainder Theorem
         * - Maps each input index to optimal (row, column) position
         * - Minimizes twiddle factor complexity
         *
         * Reverse generation (rgen=true): Use direct row-major mapping
         * - Simpler mapping for easier hardware implementation
         * - May have different memory access patterns
         */
        val load_tab: Array[Int] =
            if (rgen) {
                // Reverse generation: row-major mapping
                (for (s <- 0 until second.rN; f <- 0 until first.rN) yield
                    (second.load_tab(s) * first.rN + first.load_tab(f) * second.rN) % N).toArray
            } else {
                // Forward generation: Chinese Remainder Theorem mapping
                (for (s <- 0 until second.rN; f <- 0 until first.rN) yield
                    base.chinese(first.load_tab(f), second.load_tab(s))).toArray
            }

        /**
         * Store table generation for PSB decomposition
         *
         * Maps two-dimensional results back to output indices.
         * Uses opposite strategy from load table to ensure correct
         * input-output mapping.
         *
         * Forward generation (rgen=false): Use direct row-major mapping
         * Reverse generation (rgen=true): Use Chinese Remainder Theorem
         */
        val store_tab: Array[Int] =
            if (rgen) {
                // Reverse generation: Chinese Remainder Theorem mapping
                (for (f <- 0 until first.rN; s <- 0 until second.rN) yield
                    base.chinese(first.store_tab(f), second.store_tab(s))).toArray
            } else {
                // Forward generation: row-major mapping
                (for (f <- 0 until first.rN; s <- 0 until second.rN) yield
                    (second.store_tab(s) * first.rN + first.store_tab(f) * second.rN) % N).toArray
            }

        /**
         * String representation for debugging
         * @return String showing PSB decomposition and generation mode
         */
        override def toString = s"psb($rgen): $first|$second"
    }

    /**
     * Mixed Space Base FFT decomposition class
     *
     * Implements MSB (Mixed Space Base) FFT decomposition for general factor combinations.
     * Used when N = N1 * N2 for any factors N1, N2 (not necessarily coprime).
     * This is the most general-purpose decomposition method.
     *
     * Mathematical Foundation:
     * Based on the standard Cooley-Tukey FFT decomposition:
     * X[k] = Σ[n=0 to N-1] x[n] * W_N^(k*n)
     * where k = k1 + N1*k2 and n = n1 + N2*n2
     *
     * Decomposition Strategy:
     * Unlike PSB, MSB works for any factorization without requiring coprimality.
     * The algorithm follows the classic Cooley-Tukey approach:
     * 1. N2 transforms of size N1 (inner stage)
     * 2. N1 transforms of size N2 (outer stage)
     * 3. Twiddle factors between stages
     *
     * Index Mapping:
     * Maps 1D index k to 2D coordinates (k1, k2) using:
     * - k1 = k mod N1 (inner transform index)
     * - k2 = floor(k / N1) (outer transform index)
     *
     * Advantages:
     * - Works for any factorization
     * - General-purpose and widely applicable
     * - Well-understood mathematical properties
     *
     * Tradeoffs:
     * - More complex twiddle factors than PSB
     * - Less optimal memory access patterns for coprime factors
     * - Higher computational overhead for certain factorizations
     *
     * @param first First stage FFT processor (size N1)
     * @param second Second stage FFT processor (size N2)
     * @param rgen Whether to use reverse generation mode
     * @tparam T Implementation method type
     *
     * @example
     * // Create MSB for 12-point FFT (3 * 4, not coprime)
     * val stage3 = pb[Flat](3, false)
     * val stage4 = pb[Flat](4, false)
     * val msb12 = msb[Flat](stage3, stage4, false)
     *
     * // Create MSB for 64-point FFT (8 * 8, same factors)
     * val stage8a = pb[Flat](8, false)
     * val stage8b = pb[Flat](8, false)
     * val msb64 = msb[Flat](stage8a, stage8b, false)
     */
    case class msb[T <: Method](first: ffter[Method], second: ffter[Method], rgen: Boolean) extends ffter[T] {
        /** Total transform size N = N1 * N2 */
        val N = rN

        /** FFT specification - use MSB algorithm */
        def fft = msbFFT(this)

        /**
         * Load table generation for MSB decomposition
         *
         * Maps input indices to two-dimensional access pattern based on
         * standard Cooley-Tukey index mapping.
         *
         * Forward generation (rgen=false): Column-major access
         * - Organizes data for N1-point transforms
         * - Better cache locality for inner stage computation
         *
         * Reverse generation (rgen=true): Row-major access
         * - Different memory access pattern
         * - May be more suitable for certain hardware implementations
         */
        val load_tab: Array[Int] = if (rgen) {
            // Reverse generation: row-major mapping
            (for (f <- 0 until first.rN; s <- 0 until second.rN) yield
                first.load_tab(f) + second.load_tab(s) * first.rN).toArray
        } else {
            // Forward generation: column-major mapping
            (for (s <- 0 until second.rN; f <- 0 until first.rN) yield
                first.load_tab(f) * second.rN + second.load_tab(s)).toArray
        }

        /**
         * Store table generation for MSB decomposition
         *
         * Maps two-dimensional results back to output indices.
         * Uses the inverse mapping from load table to ensure correct
         * input-output correspondence.
         *
         * The modulo operation ensures indices stay within valid range [0, N-1].
         */
        val store_tab: Array[Int] = if (rgen) {
            // Reverse generation: column-major mapping
            (for (s <- 0 until second.rN; f <- 0 until first.rN) yield
                (second.store_tab(s) + first.store_tab(f) * second.rN) % N).toArray
        } else {
            // Forward generation: row-major mapping
            (for (f <- 0 until first.rN; s <- 0 until second.rN) yield
                (second.store_tab(s) * first.rN + first.store_tab(f)) % N).toArray
        }

        /**
         * String representation for debugging
         * @return String showing MSB decomposition and generation mode
         */
        override def toString = s"msb($rgen): $first|$second"
    }

    /**
     * Winograd Prime Base FFT decomposition class
     *
     * Implements WPBP (Winograd Prime Base) FFT decomposition for prime-sized FFTs
     * using Winograd's minimal multiplication algorithms. This is the most
     * computationally efficient method for prime-sized transforms.
     *
     * Mathematical Foundation:
     * Based on Winograd's work on minimal multiplication algorithms for DFT.
     * The key insight is that prime-length DFTs can be computed using
     * convolution with fewer multiplications than the direct DFT formula.
     *
     * Algorithm Principles:
     * 1. Convert DFT to cyclic convolution using number-theoretic properties
     * 2. Apply Winograd's minimal multiplication convolution algorithms
     * 3. Use circular permutation to optimize computational structure
     * 4. Exploit symmetries in prime-length transforms
     *
     * Circular Permutation Strategy:
     * The algorithm uses a circular ordering of input/output indices to
     * optimize the computational structure and minimize the number of
     * required multiplications.
     *
     * Performance Characteristics:
     * - Minimal multiplication count for prime sizes
     * - O(N log N) complexity vs O(N²) for direct DFT
     * - Particularly effective for medium-sized primes (5, 7, 11, 13, etc.)
     * - Hardware-friendly structure with regular data flow
     *
     * @param N Transform size (must be prime)
     * @param first First stage decomposition (typically N-1 point)
     * @param second Second stage decomposition (typically small)
     * @param rgen Whether to use reverse generation mode
     * @tparam T Implementation method type
     *
     * @throws AssertionError when N is not prime
     *
     * @example
     * // Create Winograd FFT for 7-point transform
     * val sub1 = FFTGen.gen[Flat](6)  // N-1 decomposition
     * val sub2 = fft1                // Base case
     * val winograd7 = wpbp[Flat](7, sub1, sub2, false)
     */
    case class wpbp[T <: Method](N: Int, first: ffter[Method], second: ffter[Method], rgen: Boolean) extends ffter[T] {
        // Validate that N is prime
        assert(Prime.isPrime(N), s"wpbp need $N is prime")

        /** Circular ordering table for Winograd optimization */
        val oc = circularOrder(N)

        /**
         * Transform size for Winograd decomposition
         * @return N (prime transform size)
         */
        override def rN = N

        /** FFT specification - use Winograd Prime Base algorithm */
        def fft = wpbFFT(N, this)

        /**
         * String representation for debugging
         * @return String showing Winograd FFT configuration
         */
        override def toString = s"winograd prime fft($rgen) $N "

        /**
         * Load table generation for Winograd decomposition
         *
         * Implements the circular permutation strategy that optimizes
         * the Winograd algorithm structure. The key insight is that
         * the DC component (index 0) is handled separately, and the
         * remaining N-1 components are processed using circular ordering.
         *
         * Algorithm:
         * - Index 0 maps to 0 (DC component)
         * - Other indices use circular permutation: oc.tab(original_index + 1)
         * - This ordering optimizes the convolution structure
         */
        val load_tab: Array[Int] = (for (f <- 0 until N) yield
            if (f == 0) 0
            else oc.tab(first.load_tab(f - 1) + 1)
        ).toArray

        /**
         * Store table generation for Winograd decomposition
         *
         * Applies the inverse circular permutation to reconstruct the
         * proper output ordering. The DC component remains at index 0,
         * while other components are reordered back to natural order.
         *
         * Algorithm:
         * - Index 0 maps to 0 (DC component)
         * - Other indices use inverse circular permutation
         * - Restores natural frequency domain ordering
         */
        val store_tab: Array[Int] = (for (s <- 0 until N) yield
            if (s == 0) 0
            else oc.tab(second.store_tab(s - 1) + 1)
        ).toArray
    }

    /**
     * Global algorithm selection flag
     *
     * Controls whether to use Winograd algorithms for prime factor FFTs
     * or standard flattened DFT implementations. This affects the
     * decomposition strategy throughout the FFT generation process.
     *
     * Winograd Mode (wMethod = true):
     * - Use Winograd prime factor algorithms when possible
     * - Minimal multiplication count for prime sizes
     * - More complex implementation but better performance
     *
     * Flatten DFT Mode (wMethod = false):
     * - Use standard DFT decomposition
     * - Simpler implementation with more multiplications
     * - Better for small prime sizes or when simplicity is preferred
     */
    var wMethod = false

    /**
     * Enable Winograd algorithm mode
     *
     * Sets the global flag to use Winograd algorithms for prime factor
     * FFTs. This affects all subsequent FFT generation calls.
     *
     * Effect:
     * - Prime factors will use WPB (Winograd Prime Base) when N > 2
     * - Reduces multiplication count significantly for medium primes
     * - Increases implementation complexity
     *
     * Usage:
     * Call this method before generating FFT structures if performance
     * is critical and implementation complexity is acceptable.
     *
     * @example
     * FFTGen.Winograd
     * val fft = FFTGen.gen[Flat](35)  // Will use Winograd for factor 7
     */
    def Winograd: Unit = { wMethod = true }

    /**
     * Enable Flatten DFT mode
     *
     * Sets the global flag to use standard flattened DFT algorithms
     * for all prime factors. This is the default mode.
     *
     * Effect:
     * - Prime factors will use PB (Prime Base) direct DFT
     * - Simpler implementation with more multiplications
     * - Better for small primes or when implementation simplicity is preferred
     *
     * Usage:
     * Call this method before generating FFT structures if simplicity
     * is preferred over optimal performance.
     *
     * @example
     * FFTGen.FlattenDFT
     * val fft = FFTGen.gen[Flat](35)  // Will use direct DFT for factor 7
     */
    def FlattenDFT: Unit = { wMethod = false }

    
    /**
     * Choose optimal factor for FFT decomposition
     *
     * Selects the best factorization strategy based on the input size N.
     * Uses prime factorization to determine optimal decomposition approach
     * that minimizes intermediate result sizes and computational complexity.
     *
     * Algorithm Strategy:
     * 1. Factor N into prime components using Prime.repesent()
     * 2. For single prime factor, use power-based decomposition
     * 3. For multiple factors, select optimal combination using heuristic
     *
     * Single Prime Factor Strategy:
     * - If N = 2^k with k > 2: use 4^(k/4) as factor (maps well to hardware)
     * - Otherwise: use prime^(exponent/2) as factor (balanced decomposition)
     * - This creates roughly square decompositions for efficiency
     *
     * Multiple Factor Strategy:
     * - Sort prime factors by size (largest first)
     * - Use recursive heuristic to find optimal factor combination
     * - Goal: minimize intermediate result sizes to prevent overflow
     * - Constraint: keep partial products below sqrt(N) when possible
     *
     * Selection Criteria:
     * - Minimize intermediate result sizes
     * - Balance computational load between stages
     * - Prefer factors that lead to efficient hardware implementation
     * - Avoid factors that create excessive intermediate values
     *
     * @param n The FFT size to factorize
     * @return Selected factor for decomposition (1 <= result <= n)
     *
     * @throws IllegalArgumentException When n cannot be factorized effectively
     *
     * @complexity O(log n) for prime factorization + O(factors) for heuristic search
     *
     * @example
     * val factor64 = FFTGen.choose_factor(64)  // Returns 8 (64 = 8 * 8)
     * val factor30 = FFTGen.choose_factor(30)  // Returns 6 (30 = 6 * 5)
     * val factor49 = FFTGen.choose_factor(49)  // Returns 7 (49 = 7 * 7)
     */
    def choose_factor(n: Int): Int = {
        // Get prime factorization: List[(prime, exponent)]
        val rn = Prime.repesent(n).toList

        if (rn.size == 1) {
            // Single prime factor case
            val (prime, exponent) = rn(0)

            if (prime == 2 && exponent % 2 == 0 && exponent > 2) {
                // Power of 2 case: use 4^(k/4) for better hardware mapping
                // Example: 2^8 = 256 -> 4^2 = 16 (256 = 16 * 16)
                FFTMisc.pow(4, exponent / 4).toInt
            } else {
                // General prime case: use prime^(exponent/2)
                // Example: 3^6 = 729 -> 3^3 = 27 (729 = 27 * 27)
                FFTMisc.pow(prime, exponent / 2).toInt
            }
        } else {
            // Multiple factors case: find optimal combination
            // Convert to actual factor values and sort by size (largest first)
            val np = rn.map { case (k, v) => FFTMisc.pow(k, v) }.sortWith(_ > _)

            /**
             * Recursive heuristic to find optimal factor combination
             *
             * This function tries to build up a factor by multiplying primes
             * while keeping the intermediate result below sqrt(n).
             *
             * @param a Current accumulated product
             * @param b Remaining factors to consider
             * @return Optimal factor value
             */
            def check(a: BigInt, b: List[BigInt]): BigInt = b match {
                case head :: tl =>
                    val c = a * head
                    tl match {
                        case Nil =>
                            // Last factor: check if product squared is less than n
                            if ((c * c) < BigInt(n)) c else a
                        case _ =>
                            // More factors: continue recursion if product is small enough
                            if ((c * c) < BigInt(n)) check(c, tl) else check(a, tl.tail)
                    }
                case Nil => a
            }

            // Start with largest factor and apply heuristic
            check(np.head, np.tail).toInt
        }
    }

    /**
     * Generate forward FFT decomposition structure
     *
     * Creates an optimal decomposition tree for computing forward FFT of size N.
     * The method automatically selects the best decomposition strategy based on
     * the input size and current algorithm preferences.
     *
     * Decomposition Strategy Selection:
     * 1. Prime Factor Case: If N is prime, use appropriate prime-size algorithm
     * 2. Composite Case: Factor N = i * o and apply Cooley-Tukey decomposition
     * 3. Algorithm Selection: Choose PSB or MSB based on factor coprimality
     * 4. Prime Handling: Use Winograd or standard DFT based on wMethod setting
     *
     * Algorithm Flow:
     * 1. Factor selection: Choose i (inner factor) and o (outer factor)
     * 2. Prime handling: Use specialized algorithms for prime sizes
     * 3. Recursive decomposition: Break down factors until base cases
     * 4. Method combination: Combine sub-transforms using PSB or MSB
     *
     * Prime Factor Handling:
     * - For N = 2: Use direct DFT (always optimal)
     * - For N > 2 prime: Use Winograd if wMethod=true, else direct DFT
     * - For N = prime + 1: Use recursive Winograd construction
     *
     * @param N FFT transform size to generate
     * @tparam T Implementation method type (Flat, Mixed, Streamed)
     * @return FFT decomposition structure for computing the transform
     *
     * @throws IllegalArgumentException When N cannot be decomposed effectively
     *
     * @complexity O(log N) for factorization + recursion depth for tree generation
     *
     * @example
     * // Generate 64-point FFT decomposition
     * val fft64 = FFTGen.gen[Flat](64)
     * // Result: psb(psb(pb(8), pb(8)), psb(pb(8), pb(8))) for 64 = 8*8*8
     *
     * // Generate 30-point FFT decomposition
     * val fft30 = FFTGen.gen[Mixed](30)
     * // Result: psb(pb(5), pb(6)) for 30 = 5*6 (coprime factors)
     */
    def gen[T <: Method](N: Int): ffter[T] = {
        /**
         * Factor selection and initial decomposition
         *
         * If N is prime, use N as factor (will be handled by prime algorithm)
         * Otherwise, use choose_factor() to find optimal decomposition factor
         */
        val i = if (Prime.isPrime(N)) N else choose_factor(N)
        val o = N / i

        /**
         * Prime factor handling function
         *
         * Handles the base case when a factor is prime or when the size is small.
         * Uses different strategies based on size and algorithm preferences.
         *
         * @param n Size to handle as prime factor
         * @return FFT structure for prime-size transform
         */
        def pp(n: Int) = {
            if (n == 2 || !wMethod) {
                // Small size or Winograd disabled: use basic prime base DFT
                pb[T](n, false)
            } else {
                // Winograd enabled for larger prime sizes: use recursive construction
                val first = gen[T](n - 1).asInstanceOf[ffter[Method]]
                val second = rgen[T](n - 1).asInstanceOf[ffter[Method]]
                wpbp[T](n, first, second, false)
            }
        }

        // Base case handling: when one factor is 1, handle entire transform as prime
        if (i == 1) pp(N)
        else if (o == 1) pp(N)

        // Composite case: apply Cooley-Tukey decomposition
        else {
            // Generate decomposition for inner factor
            val ib = gen[T](i).asInstanceOf[ffter[Method]]
            // Generate decomposition for outer factor
            val ob = gen[T](o).asInstanceOf[ffter[Method]]

            /**
             * Choose combination method based on factor properties
             *
             * - PSB (Prime Space Base): Used when factors are coprime (gcd = 1)
             *   More efficient memory access patterns
             * - MSB (Mixed Space Base): Used when factors share common divisors
             *   More general case, slightly less efficient
             */
            if (Prime.gcd(i, o) == 1) psb(ib, ob, false)  // Coprime factors: use PSB
            else msb(ib, ob, false)                       // Shared factors: use MSB
        }
    }

    /**
     * Generate reverse FFT decomposition structure
     *
     * Creates an optimal decomposition tree for computing reverse FFT (IFFT) of size N.
     * Similar to gen() but uses reverse generation mode which affects twiddle factor
     * signs and data ordering patterns.
     *
     * Key Differences from gen():
     * - Uses reverse generation mode (rgen=true) in all subcomponents
     * - Affects twiddle factor signs (conjugate for inverse FFT)
     * - Different data permutation patterns for optimal IFFT computation
     * - Same factorization strategy as forward FFT
     *
     * Reverse Generation Effects:
     * - Twiddle factors use opposite signs: W_N^(-k*n) instead of W_N^(k*n)
     * - Load/store tables use different permutations
     * - Final output scaling may be required (1/N factor)
     *
     * @param N FFT transform size to generate
     * @tparam T Implementation method type (Flat, Mixed, Streamed)
     * @return FFT decomposition structure for computing inverse transform
     *
     * @see gen for forward FFT generation
     *
     * @example
     * // Generate 64-point inverse FFT decomposition
     * val ifft64 = FFTGen.rgen[Flat](64)
     * // Similar structure to gen() but with reverse generation parameters
     */
    def rgen[T <: Method](N: Int): ffter[T] = {
        /**
         * Factor selection for reverse generation
         *
         * Note: Uses same factor selection as forward FFT but with
         * different ordering strategy (i and o are swapped compared to gen)
         */
        val o = if (Prime.isPrime(N)) N else choose_factor(N)
        val i = N / o

        /**
         * Prime factor handling for reverse generation
         *
         * Same logic as pp() in gen() but with reverse generation flag
         */
        def pp(n: Int) = {
            if (n == 2 || !wMethod) {
                // Small size or Winograd disabled: use basic prime base DFT
                pb[T](n, true)
            } else {
                // Winograd enabled for larger prime sizes: use recursive construction
                // Note: Mixes gen() and rgen() for optimal structure
                val first = gen[T](n - 1).asInstanceOf[ffter[Method]]
                val second = rgen[T](n - 1).asInstanceOf[ffter[Method]]
                wpbp[T](n, first, second, true)
            }
        }

        // Base case handling: same logic as gen() but with reverse generation
        if (i == 1) pp(N)
        else if (o == 1) pp(N)

        // Composite case: apply Cooley-Tukey decomposition with reverse generation
        else {
            // Generate decomposition with reverse generation mode
            val ib = rgen[T](i).asInstanceOf[ffter[Method]]
            val ob = rgen[T](o).asInstanceOf[ffter[Method]]

            // Same combination method selection as gen() but with reverse generation
            if (Prime.gcd(i, o) == 1) psb(ib, ob, true)  // Coprime factors: use PSB
            else msb(ib, ob, true)                       // Shared factors: use MSB
        }
    }
}

